module sw_case6_mod

  use const_mod

  implicit none

  private

  public sw_case6

  real(8), parameter :: radius  = 6371.220d3
  real(8), parameter :: gravity = 9.80616d0
  real(8), parameter :: omega   = 7.292d-5
  real(8), parameter :: R       = 4.0d0
  real(8), parameter :: omg     = 7.848d-6
  real(8), parameter :: gh0     = 8.0d3 * gravity

contains

  subroutine sw_case6(lon, lat, u, v, gh, ghs, g, f, radi)

    real(8), intent(in) :: lon
    real(8), intent(in) :: lat
    real(8), intent(out), optional :: u
    real(8), intent(out), optional :: v
    real(8), intent(out), optional :: gh
    real(8), intent(out), optional :: ghs
    real(8), intent(out), optional :: g
    real(8), intent(out), optional :: f
    real(8), intent(out), optional :: radi

    real(8) cos_lon, sin_lon, cos_lat, sin_lat
    real(8) a, b, c

    if(present(radi))radi = radius
    
    cos_lon = cos(lon)
    sin_lon = sin(lon)
    cos_lat = cos(lat)
    sin_lat = sin(lat)

    if (present(u)) then
      a = cos_lat
      b = R * cos_lat**(R - 1) * sin_lat**2 * cos(R * lon)
      c = cos_lat**(R + 1) * cos(R * lon)
      u = radius * omg * (a + b - c)
    end if
    if (present(v)) then
      a = R * cos_lat**(R - 1) * sin_lat * sin(R * lon)
      v = -radius * omg * a
    end if
    if (present(gh)) then
      a = 0.5 * omg * (2 * omega + omg) * cos_lat**2 + &
        0.25 * omg**2 * ((R + 1) * cos_lat**(2 * R + 2) + (2 * R**2 - R - 2) * cos_lat**(2 * R) - 2 * R**2 * cos_lat**(2 * R - 2))
      b = 2 * (omega + omg) * omg * cos_lat**R * &
        (R**2 + 2 * R + 2 - (R + 1)**2 * cos_lat**2) / (R + 1) / (R + 2)
      c = 0.25 * omg**2 * cos_lat**(2 * R) * ((R + 1) * cos_lat**2 - R - 2)
      gh = gh0 + radius**2 * (a + b * cos(R * lon) + c * cos(2 * R *lon))
    end if
    if (present(ghs)) ghs = 0
    if (present(g  )) g   = gravity
    if (present(f  )) f   = 2 * omega * sin_lat

  end subroutine sw_case6

end module sw_case6_mod
